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Mimicking  celestial  mechanics  in  metamaterials 

Dentcho  A.  Genov1'2,  Shuang  Zhang1  and  Xiang  Zhang1,3* * 

Einstein's  general  theory  of  relativity  establishes  equality  between  matter-energy  density  and  the  curvature  of  spacetime.  As 
a  result,  light  and  matter  follow  natural  paths  in  the  inherent  spacetime  and  may  experience  bending  and  trapping  in  a  specific 
region  of  space.  So  far,  the  interaction  of  light  and  matter  with  curved  spacetime  has  been  predominantly  studied  theoretically 
and  through  astronomical  observations.  Here,  we  propose  to  link  the  newly  emerged  field  of  artificial  optical  materials  to  that  of 
celestial  mechanics,  thus  opening  the  way  to  investigate  light  phenomena  reminiscent  of  orbital  motion,  strange  attractors  and 
chaos,  in  a  controlled  laboratory  environment.  The  optical-mechanical  analogy  enables  direct  studies  of  critical  light/matter 
behaviour  around  massive  celestial  bodies  and,  on  the  other  hand,  points  towards  the  design  of  novel  optical  cavities  and  photon 
traps  for  application  in  microscopic  devices  and  lasers  systems. 


The  possibility  to  precisely  control  the  flow  of  light  by  designing 
the  microscopic  response  of  an  artificial  optical  material  has 
attracted  great  interest  in  the  field  of  optics.  This  ongoing 
revolution,  facilitated  by  the  advances  in  nanotechnology,  has 
enabled  the  manifestation  of  exciting  effects  such  as  negative 
refraction1,2,  electromagnetic  invisibility  devices3-5  and  microscopy 
with  super-resolution1,6.  The  basis  for  some  of  these  phenomena 
is  the  equivalence  between  light  propagation  in  curved  space  and 
in  a  locally  engineered  optical  material.  Interestingly,  a  similar 
behaviour  exists  in  the  general  theory  of  relativity,  where  the 
presence  of  matter-energy  densities  results  in  curved  spacetime 
and  complex  motion  of  both  matter  and  light7,8.  In  the  classical 
interpretation,  this  fundamental  behaviour  is  known  as  the  optical- 
mechanical  analogy  and  is  revealed  through  the  least- action 
principles  in  mechanics,  determining  how  a  particle  moves  in  an 
arbitrary  potential9,  and  the  Fermat  principle  in  optics,  describing 
the  ray  propagation  in  an  inhomogeneous  media10. 

The  optical-mechanical  analogy  provides  a  rather  useful  link 
between  the  study  of  light  propagation  in  inhomogeneous  media 
and  the  motion  of  massive  bodies  or  light  in  gravitational  potentials. 
Surprisingly,  a  direct  mapping  of  the  celestial  phenomena  by  ob¬ 
serving  photon  motion  in  a  controlled  laboratory  environment  has 
not  been  studied  so  far  in  an  actual  experiment.  Here,  we  propose 
realistic  optical  materials  that  facilitate  periodic,  quasi-periodic 
and  chaotic  light  motion  inherent  to  celestial  objects  subjected 
to  complex  gravitational  fields.  This  dense  optical  media  (DOM) 
could  be  readily  achieved  with  the  current  technology,  within 
the  framework  of  artificial  optical  metamaterials.  Furthermore, 
we  introduce  a  new  class  of  specifically  designed  DOM  in  the 
form  of  continuous -index  photon  traps  (CIPTs)  that  can  serve 
as  novel  broadband,  radiation-free  and  thus  ‘perfect5  cavities.  The 
dynamics  of  the  electromagnetic  energy  confinement  in  the  CIPTs 
is  demonstrated  through  full-wave  calculations  in  the  linear  and 
nonlinear  regimes.  Possible  mechanisms  for  coupling  light  into  the 
CIPTs  and  specific  composite  media  are  discussed. 

Optical  attractors  and  photonic  black  holes 

One  of  the  most  fascinating  predictions  of  the  general  theory  of 
relativity  is  the  bending  of  light  that  passes  near  massive  celestial 
objects  such  as  stars,  nebulas  or  galaxies.  This  effect  constituted 


one  of  the  first  pieces  of  evidence  for  the  validity  of  Einstein’s 
theory  and  provided  a  glimpse  into  the  interesting  predictions 
that  were  to  follow.  As  shown,  first  by  Schwarzschild  and  then  by 
others,  when  the  mass  densities  of  celestial  bodies  reach  a  certain 
critical  value  (for  example,  through  gravitation  collapse),  objects 
termed  black  holes  form  in  space.  Those  entities  fall  under  a  more 
general  class  of  dynamic  systems,  where  a  particular  point,  curve 
or  manifold  in  space  acts  as  an  attractor  of  both  matter  and  light. 
Such  systems  are  of  great  interest  for  science  not  only  in  terms  of 
the  fundamental  studies  of  light-matter  interactions  but  also  for 
possible  applications  in  optical  devices  that  control,  slow  and  trap 
light.  It  is  thus  important  to  investigate  the  associated  phenomena 
under  controlled  laboratory  conditions. 

In  the  general  theory  of  relativity,  a  complex  particle/photon 
motion  is  observed  whenever  the  inherent  spacetime  is  described  by 
a  metric  gy  that  depends  on  the  spatial  coordinates  x  =  {x1,^2,^3} 
and  universal  time  t.  In  such  curved,  non-static  spacetimes,  the 
propagation  of  matter  and  light  rays  follows  the  natural  geodesic 
lines  and  is  described  by  the  Lagrangian 

1  ,1 

L=-g00(x,t)i  - -gij(x,t)xlxJ  (1) 

where  the  derivatives  are  taken  over  an  arbitrary  affine  parameter 
and  natural  units  have  been  adopted.  Here,  we  investigate  the 
prospect  to  mimic  the  effect  of  curved  spacetimes  (equation  ( 1))  on 
matter  and  light  in  the  laboratory.  To  achieve  this,  we  recognize  that 
complex  light  motion  is  also  observed  in  metamaterials  described 
by  inhomogeneous  permittivities  and  permeabilities  and  could  be 
related  to  light  dynamics  in  curved  space  through  the  invariance 
of  Maxwell’s  equations  under  coordinate  transformations4,5,11. 
Metamaterials  exhibiting  complex  electric  and  magnetic  responses 
have  been  the  focus  of  recent  efforts  to  create  negative- refractive - 
index  media1,2,  electromagnetic  cloaking2-4  and  as  concentrators 
of  light12.  The  above  systems,  however,  suffer  from  high  intrinsic 
loss  and  a  narrow  frequency  range  of  operation  and  thus  cannot 
be  considered  as  prospective  media  to  simulate  light  motion 
in  a  curved  spacetime  vacuum.  Although  the  introduction  of 
new  metamaterials13  may  improve  the  overall  performance,  it  is 
important  to  rely  only  on  optical  materials  that  are  non-dissipative 
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Figure  1 1  Optical  attractors  and  PBHs.  a,  Curved  space  corresponding  to  a  PBH  with  goo  =  1  showing  regions  of  spatial  stretching  and  compression. 
Depending  on  their  initial  angular  momentum  (impact  parameter  s),  the  rays  approaching  the  PBH  will  be  deflected  (s  >  b ),  attracted  (s  <  b)  or  settle 
(s  =  b),  on  an  unstable  orbit.  b,c,  A  set  of  such  critical  photon  trajectories  that  fall  into  the  singularity  for  a  =  0  (b)  or  approach  the  unstable  orbit  at  a  =  1 
(c)  for  an  impact  parameter  b  =  2.  d,e,  The  FDFD  calculation  of  the  magnetic  field  density  (the  incident  magnetic  field  is  |Ho|  =  1),  corresponding  to  the 
critical  trajectories  with  b  =  1.5a  =  15  gm  (d)  and  b  =  0.5 a  =  5  gm  (e).  In  the  calculations,  very  high  energy  densities  are  obtained  close  to  the  central 
attractor  with  the  maximum  value  restricted  only  by  the  finite  length  scale  that  could  be  numerically  investigated.  The  incident  wavelength  is 
k  =  0.5  pirn  and  the  Gaussian  beam  width  is  w  =  4k. 


and  non-dispersive  in  general.  For  that,  we  consider  a  class  of 
centrally  symmetric  metrics  that  can  be  written  under  coordinate 
transformation  in  an  isotropic  form  (see  the  Methods  section).  For 
such  curved  spacetimes,  light  behaves  in  space  similarly  as  in  an 
optical  media  with  an  effective  refraction  index 

n  =  (g/goo)1/2  (2) 

where  g  =  gn  —  gn  —  £33-  The  equivalence  manifested  in 
equation  (2)  is  the  basis  for  the  studies  of  celestial  phenomena  in  an 
actual  table-top  experiment.  It  enables  the  use  of  non-dissipative 
and  non-dispersive  dielectric  materials.  To  demonstrate  this,  we 
propose  a  two -parameter  family  of  isotropic  centrally  symmetric 
metrics  g  —  goo((b/r)2  +  (1  —  a/r)2),  where  r  =  |x|  is  the  radial 
coordinate  and  a  and  b  are  constants;  constraint  lim^oogoo  =  1 
is  imposed  to  recover  the  flat  space  at  large  distances.  The  choice 
of  this  particular  spacetime,  which  is  schematically  represented  in 
Fig.  la,  has  a  number  of  important  rationales.  It  exhibits  a  physical 
singularity  at  r  =  0  and  as  shown  below  unstable  photon  orbits 
are  supported  at  r  —  a,  defining  a  photon  sphere  for  the  system14. 
These  characteristics  are  reminiscent  of  a  gravitational  black  hole 
in  the  general  theory  of  relativity  although  the  proposed  metrics 
are  not  solutions  to  the  Einstein  field  equations  in  a  vacuum. 
Notwithstanding,  here  we  refer  to  the  considered  system  as  a 
photonic  black  hole  (PBH),  or  a  spacetime  where  a  particular  point 
in  space  acts  as  an  attractor  of  electromagnetic  radiation. 

Figure  lb,c  demonstrates  photon  trajectories  incident  on  the 
PBH  with  an  impact  parameter  s,  representing  the  minimal  distance 


to  the  centre  provided  the  photon  is  not  deflected,  that  is  equal 
to  a  critical  value  b.  The  Lagrangian  formalism  based  on  the 
curved  spacetime  equation  (1)  predicts  that  the  impinging  rays  will 
inevitably  reach  their  destiny  at  the  central  attractor/singularity 
(Fig.  lb),  or  asymptotically  approach  the  photon  sphere  at  r  —  a 
(Fig.  lc).  The  spiral  trajectories  represent  a  critical  point  in  the 
PBH  phase  space  (see  Supplementary  Information)  such  that  all 
photons  approaching  with  impact  factors  less  than  or  equal  to 
the  critical  parameter  b  will  never  escape  the  system.  Figure  ld,e 
shows  full-wave  simulations  of  Gaussian  beams  incident  on  a 
DOM  with  a  refraction  index  profile  described  by  equation  (2) 
that  maps  the  PBH  space.  The  incident  beams  have  propagation 
directions  and  initial  positions  set  on  the  PBH  critical  trajectories 
(as  in  Fig.  lb,c),  and  the  light  is  concurrently  torn  apart,  with 
half  the  beams  being  scattered  into  infinity  following  hyperbolical 
paths,  whereas  the  rest  spirals  into  the  central  attractor.  This  ray 
behaviour  is  consistent  with  the  expected  ray  trajectories  in  the 
original  curved  spacetime,  but  demonstrated  in  the  equivalent 
DOM.  In  addition  to  the  predicted  ray  behaviour,  we  observe 
electromagnetic  energy  density  fluctuations  (jets)  leaving  the  system 
in  radial  directions  (see  Fig.  le).  This  effect  is  due  to  the  wave  nature 
of  light  interaction  with  the  inhomogeneous  refraction  index,  and 
affects  the  electromagnetic  rays  that  approach  the  PBH  with  impact 
factors  slightly  above  the  critical  value.  At  the  centre  of  the  system, 
owing  to  the  spacial  singularity  inherent  to  the  PBH  metrics, 
the  effective  refractive  index  approaches  infinity.  The  singularly 
corresponds  to  an  increasing  stretching  of  space  (see  Fig.  la),  which 
acts  as  an  effective  ‘potential  well’  attracting  light.  Concurrently,  the 
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optical  attractor  can  be  considered  as  an  effective  ‘cavity5,  decreasing 
both  the  group  velocity  and  wavelength,  which  in  turns  gives  rise  to 
exceptionally  high  local  energy  densities. 

Although  the  singularity  in  the  refractive  index  profile  of  the 
ideal  PBH  is  of  fundamental  interest  for  science,  it  is  probably 
not  feasible  in  practice.  Such  index  profiles  or  the  associated  low 
group  velocities  have  been  discussed  in  the  literature  and  specific 
systems  have  been  suggested  such  as  high-velocity  vortex  flow15,16, 
microstructured  nonlinear  fibres17,  metamaterial  waveguides18, 
metal  nanoparticles19  and  Bose-Einstein  condensates20.  Further¬ 
more,  close  to  the  centre  singularity,  the  ray  picture  fails  and  one 
cannot  trust  the  predictions  as  per  Fig.  lb,c.  However,  experimental 
observation  of  the  main  optical  phenomena  related  to  systems 
such  as  PBH  or  gravitational  black  holes  does  not  necessitate 
exceptionally  large  refraction  indices,  or  extremely  low  phase/ group 
velocities.  In  fact,  next  we  show  that  the  underlying  electromagnetic 
phenomena  can  be  demonstrated  using  DOM  in  the  form  of  simple 
binary  metal-dielectric  or  pure  dielectric  composites. 

Figure  2  shows  finite  refractive  index  profiles  n  mapped  to 
that  of  a  PBH,  obtained  by  varying  the  volume  fraction  p  of 
copper  metal  nanoparticles  suspended  in  air  (Fig.  2a)  or  air  holes 
in  a  GaAs  substrate  (Fig.  2b).  These  composite  materials  can  be 
nano-engineered  through  controlled  deposition  of  metal  particles 
in  a  dielectric  host  or  growth  of  deep  subwavelength  in  size  air 
pockets  in  high-index  materials.  The  electromagnetic  phenomena 
associated  with  these  realistic  systems  are  shown  in  Fig.  2c, d, 
respectively.  Despite  the  finite  indices,  the  air-copper  mixture 
can  provide  a  complete  picture  of  the  far-field  scattering  profile 
inherent  to  a  PBH  (Fig.  2c),  and  photon  ‘trapping5  into  the  unstable 
orbit  (at  the  photon  sphere)  as  well  as  critical  behaviour  at  lower 
impact  factors  can  be  demonstrated  in  the  air-GaAs  composite 
(Fig.  2d).  Specifically,  at  sufficiently  large  distances  from  the  PBH 
(r  max(a,  b )),  the  light  rays  are  deflected  by  an  angle  </>  =  —2 a/s 
(Fig.  2c)  that  has  the  same  dependence  on  the  impact  parameter  s  as 
that  of  the  Einstein  angle  for  gravitational  lensing,  but  with  reversed 
sign8.  This  results  in  a  diverging  energy  flux  at  large  distances  with 
maximum  field  intensity  observed  at  ym  =  (8a|x|)1/2.  As  metals 
are  highly  dispersive,  the  air-copper  system  will  have  a  narrow 
operational  bandwidth  around  the  plasma  resonance  wavelength 
A.p  =  0.27  pm,  and  the  incident  light  will  suffer  energy  loss,  which 
although  much  lower  compared  with  current  metamaterials  is  still 
considerable  (the  imaginary  part  of  the  effective  index  is  shown  in 
Fig.  2a).  To  overcome  these  obstacles,  one  can  instead  design  a  pure 
dielectric  media  based  on  air-GaAs  composites.  Owing  to  the  larger 
refractive  index  of  the  semiconductor  (nmax  =  nGa As  =  3.25;  ref.  21), 
for  this  system,  it  is  even  possible  to  follow  the  light  propagation 
in  the  region  of  space  below  the  photon  sphere  with  significantly 
reduced  intrinsic  losses  and  for  a  wide  range  of  frequencies  (below 
the  GaAs  bandgap).  This  is  shown  in  Fig.  2d,  where  four  Gaussian 
beams  are  sent  into  the  system  with  initial  angular  moments  and 
positions  centred  on  critical  trajectories.  All  beams  asymptotically 
approach  the  photon  sphere  where  they  are  ripped  apart,  with  half 
the  energy  being  scattered  into  the  far-field,  and  the  rest  collapses 
towards  the  centre.  The  electromagnetic  density  fluctuations  or 
‘jets5  of  energy  leaving  the  photon  sphere  are  clearly  seen.  To 
guarantee  that  this  new  phenomenon  is  not  due  to  the  truncated 
region  below  the  photon  sphere,  an  index-matched  absorption  core 
is  used.  The  absorption  core  has  the  role  of  the  central  attractor 
and  assures  that  all  scattering  and  wave  phenomena  are  properly 
represented  by  the  DOM. 

Overall,  the  DOM  as  per  Fig.  2c  provides  the  means  to  study 
electromagnetic  wave  interactions  with  critical  points  in  spacetime 
(in  particular  the  photon  sphere)  in  a  relatively  simple  table- 
top  experiment.  Such  critical  (saddle)  points  have  an  important 
role  in  establishing  the  overall  system  dynamics  and  can  serve 
as  a  source  of  instability  and  chaos.  The  full  range  of  DOM’s 
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Figure  2  |  Mimicking  the  PBH  electromagnetic  phenomenon  in  the 
laboratory,  a-d,  Mixtures  of  low-concentration  copper  metal  nanoparticles 
in  air  (a)  and  air-GaAs  composite  media  (b)  can  be  used  to  provide  the 
far-field  scattering  of  the  PBH  with  a  scattering  angle  0  =  tan_1(ym/|x|)  (c) 
and  critical  behaviour  in  the  vicinity  of  the  photon  sphere  (d),  respectively. 
The  effective  refractive  indices  of  the  composites  are  presented  with  solid 
and  short-dashed  lines  for  the  real  and  imaginary  parts,  respectively,  and 
long-dashed  lines  are  used  for  the  metal/air  volume  fractions.  The  shaded 
areas  in  a  and  b  represent  the  region  of  space  where  the  refraction  indices 
are  truncated.  In  the  FDFD  calculations,  we  have  set  c,  o  =  0.25  pm, 
b  =  1.5a  and  d,  a  =  30  pm,  b  =  50  pm.  The  incident  wavelengths  and  beam 
widths  are  A.  =  0.24  pm  (w  =  20A)  and  A  =  1.55  pm  (w  =  8A),  respectively. 
For  the  air-GaAs  composite,  the  absorption  coefficient  is  fixed  at  5  cm-1 
(the  intrinsic  semiconductor  absorption  for  frequency  below  the 
bandgap)21  everywhere  in  space  except  within  the  truncation  region  where 
it  is  increased  to  100  cm-1  to  simulate  the  energy  sink  represented  by  the 
central  attractor. 


applications  including  concentrating  and  storing  electromagnetic 
energy  are  discussed  next. 

Photon  traps,  orbital  motion  and  high-Q-factor  cavities 

The  equivalence  between  light  ray  motion  around  massive  celestial 
objects  and  in  an  artificially  engineered  optical  medium  has 
strong  implications  as  the  general  formalism,  represented  by  the 
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Figure  3  |  CIPTs  based  on  air-GalnAsP  composite  media.  a,b  Refractive  index  profile  (a)  and  volume  fraction  (b)  of  the  GalnAsP.  c,  FDFD  simulations  of 
bound  and  stable  light  propagation  in  the  form  of  A,  a  circular  orbit  with  period  2jt  and  B,  a  'Rosetta'  type  of  orbit  with  period  4jt.  In  the  calculations,  a 
linear  source  with  width  w  =  6X  and  wavelength  X  =  1.55  pm  is  used,  d-f,  Poincare  maps  (Pj)  corresponding  to  all  possible  photon  trajectories  for  linear  (d) 
and  nonlinear  (e,f)  refractive  indices.  In  the  linear  case,  the  photon  motion  follows  stable  periodic  and  quasi-periodic  trajectories  for  all  initial  conditions. 
With  the  introduction  of  nonlinearity,  the  stability  of  the  system  is  disturbed  with  some  trajectories  passing  through  the  cavity  boundaries  and  leaving  the 
system,  f,  For  above-critical  values  of  the  nonlinear  index  An\^i/nmax  =  0.2,  chaotic  behaviour  is  observed  similarly  to  the  three-body  problem  in  celestial 
mechanics.  In  generating  the  Poincare  maps,  the  scale-invariant  parameters  r/a  and  r/a  are  sampled  at  times  tm  =  2nm/co^i  that  are  an  integer  number  m 
of  the  nonlinear  period  (or  T-period  map),  and  the  derivative  is  taken  with  respect  to  the  azimuth  angle  <p. 


Lagrangian  equation  (1),  also  describes  restricted  motion  of  matter 
in  potential  fields.  Namely,  under  certain  conditions  the  dynamics 
of  matter  in  complex  gravitational  potentials  can  be  replicated  by 
the  propagation  of  light  in  DOM.  Here,  we  introduce  an  interesting 
prospective,  with  significant  practical  importance,  which  is  to 
mimic  planetary  motion  in  the  form  of  a  stable  and  confined 
photon  trajectory. 

According  to  Bertrand’s  theorem22,  a  general  particle  orbit 
in  classical  mechanics  is  stable  and  closed  under  any  perturba¬ 
tions  only  in  two  types  of  central  potential:  Kepler’s  problem 
U  —  —m/r  <  0  and  the  harmonic  oscillator  U  —  mr 2  >  0,  where  m 
is  the  mass.  The  theorem,  however,  demands  further  restrictions  in 
optics  where  it  can  be  shown  that  Kepler’s  potential  no  longer  sup¬ 
ports  stable  and  closed  circular  or  elliptical  orbits  for  photons.  This 
conclusion  follows  from  the  kinematic  analogy10,  which  casts  the 
dynamic  eikonal  (ray)  equation  in  a  classical  mechanical  form  with 
effective  potential  and  kinetic  energies  U  —  —K  —  —n2/ 2  directly  re¬ 
lated  to  the  refraction  index.  Thus,  in  contrast  to  a  particle  that  may 
have  positive  or  negative  total  energy  U +K,  photons  are  equivalent 
to  ‘zero’  energy  states,  which  also  follows  from  the  null  geodetic 
condition  L  0,  or  the  finite  speed  of  light.  As  a  consequence,  only 
parabolic  solutions  are  supported  in  this  particular  case. 


Although  Bertrand’s  theorem  states  that  bound  and  exponen¬ 
tially  stable  photon  orbits  that  are  closed  under  any  perturbations 
may  not  exist  in  centrally  symmetric  media,  it  does  not  preclude 
bound  and  stable  motion  that  is  not  necessarily  closed  under 
perturbations;  that  is,  the  perturbed  trajectory  may  cover  a  finite 
volume  of  space.  Here,  we  propose  a  universal  family  of  optical 
systems  that  exhibit  those  characteristics.  The  effective  refractive 
indices  can  be  directly  obtained  from  the  stability  condition  of  the 
Lagrangian  ( 1 )  resulting  in 


where  x  (r)  is  an  arbitrary  monotonously  decreasing  function  and  A 
is  a  constant  (see  Supplementary  Information).  The  above  family  of 
refractive  indices,  which  we  refer  to  as  CIPTs,  support  stable  circular 
and  complex  ray  motion. 

To  illustrate  this  statement,  we  consider  the  most  simple 
case  where  dx/dr  =  —k/cl{k  >  0).  The  corresponding  refraction 
indices  are  then  obtained  from  equation  (3)  and  can  be  written 
as  n  —  ns(a/r)l~Ke~Kr/a ,  where  a  is  the  radius  of  the  desired 
circular  orbit  and  ns  is  an  arbitrary  scaling  factor.  Furthermore, 
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we  choose  the  k  —  2  case,  and  design  the  CIPT  as  an  air-GalnAsP 
composite.  This  specific  material  system  has  been  selected  to 
provide  operation  at  the  infrared  spectral  range,  using  the  high 
refractive  index  of  the  GalnAsP  at  low  absorption,  which  for 
frequencies  below  the  bandgap  (hcog  =  1.2  eV)  is  less  than  1  cm-1 
(ref.  21).  The  actual  index  and  air  volume  fractions  are  shown 
in  Fig.  3a, b.  Full-wave  simulations  of  two  light  trajectories  for 
linear  sources  positioned  at  rin  =  a  and  rin  =  1.62a  are  shown  in 
Fig.  3c.  According  to  the  CIPT  phase  space  shown  in  Fig.  3d,  all 
ray  trajectories  are  confined  in  space  and  the  photon  trapping  is 
achieved  without  involving  index  discontinuities  at  the  boundaries 
between  different  media.  These  bulk  type  of  cavities  where  the 
trapped  light  satisfies  the  Lyapunov  stability  condition,  constitutes  a 
substantial  improvement  compared  with  the  conventional  gradient 
refraction  index  waveguides,  and  to  the  best  of  our  knowledge,  has 
not  been  studied  in  the  literature.  The  CIPTs  may  provide  new 
directions  for  the  development  of  optical  devices  and  high-Q-factor 
cavities  for  electromagnetic  energy  storage  with  minimal  radiation 
losses.  Specifically,  coupling  CIPTs  with  a  light  source,  such  as 
a  quantum  dot,  could  result  in  enhanced  spontaneous  emission, 
strong  atom-photon  coupling  (cavity  quantum  electrodynamics), 
Raman  lasing  or  enhanced  high-order  optical  nonlinearities.  These 
phenomena  are  due  to  the  decreased  photon  decay  rate,  which 
for  CIPTs  is  r  =  C  +  rr  &  C  =  colmn/Ren ,  where  i  and  r 
stand  for  intrinsic  dissipation  and  radiation  loss,  respectively. 
The  CIPT  quality  factor  Q  =  co/T  is  thus  limited  only  by 
the  intrinsic  loss  and  can  be  on  par  or  larger  than  those  of 
photonic  crystals23  and  disc  microcavities24.  In  the  case  of  the 
air-GalnAsP  composite,  the  quality  factor  can  be  as  high  as 
106  when  taking  into  account  the  actual  material  absorption 
in  the  semiconductor21.  Furthermore,  the  Laypunovs  stability 
of  the  photon  trajectories  implies  that  finite  index  variations 
due  to  fabrication  imperfections  will  not  lead  to  photon  loss 
from  the  CIPT;  the  photons  simply  hop  to  another  close-by 
bound  trajectory  (see  Fig.  3d).  This  is  another  advantage  of  CIPTs 
in  respect  to  conventional  cavities  and  waveguides  where  light 
scattering  at  surface  or  volume  imperfections  can  lead  to  photons 
escaping  from  the  system. 

Nonlinearity  and  chaos 

The  onset  of  chaos  in  dynamic  systems  is  one  of  the  most  fascinating 
problems  in  science  and  is  observed  in  areas  as  diverse  as  molecular 
motion25,  population  dynamics26  and  optics27.  In  particular,  a 
planet  around  a  star  can  undergo  chaotic  motion  if  a  perturbation, 
such  as  another  large  planet,  is  present.  However,  owing  to  the  large 
spatial  distances  between  the  celestial  bodies,  and  the  long  periods 
involved  in  the  study  of  their  dynamics,  the  direct  observation 
of  chaotic  planetary  motion  has  been  a  challenge28.  The  use  of 
the  optical-mechanical  analogy  may  enable  such  studies  to  be 
accomplished  in  the  laboratory.  Furthermore,  in  the  case  of  the 
CIPTs,  where  trapping  of  light  occurs  with  unprecedented  lifetimes, 
the  onset  of  chaos  could  enable  access  to  those  closed  ray  trajectories 
from  the  far-field.  This  is  directly  addressed  by  involving  the  analogy 
with  the  celestial  mechanics  where  a  third-body  perturbation  can 
be  modelled  as  a  periodic  nonlinear  refraction  index  modulation; 
n(r,t)  =  hCipt(?0  +  Ahnl  sin2(&>NLf)  introduced  by  an  external 
electrical  or  electromagnetic  source. 

The  Poincare  maps  of  the  photon  motion  under  these  settings 
are  shown  in  Fig.  3e,f.  In  contrast  to  the  non-perturbed  linear  case 
(Fig.  3d),  where  all  photon  trajectories  are  closed  and  determined 
by  the  first  integral  of  the  system,  the  introduction  of  nonlinear 
harmonic  variations  brings  a  fundamentally  different  behaviour. 
For  A nNL/nmax  =  0.05  (Fig.  3e),  and  period  of  the  modulations 
T  =  2k/conl  =  6c  / a,  we  observe  substantial  expansion  of  the  origi¬ 
nal  phase  space,  the  boundary  of  which  is  marked  with  a  blue  line. 
The  original  invariant  curves  have  broken  up  into  four  separate 


domains.  As  a  result,  photons  may  pass  at  distances  r  >  rmax  =  1.74a, 
where  the  CIPT  refractive  index  truncates  to  unity.  Those  light  rays 
will  escape  on  reaching  the  outer  boundary  or  concurrently,  exter¬ 
nal  light  could  be  coupled  into  the  otherwise  perfect  cavity  (in  the 
linear  sense).  This  is  accomplished  if  the  nonlinearity  is  switched  off 
for  times  t  <  a/ c,  leading  to  a  fraction  of  the  impinging  light  being 
trapped  inside  the  system.  Similar  behaviour  could  be  observed 
even  for  moderate  nonlinearities  A nNL/nmax  <  0.01,  in  which  case  a 
smaller  fraction  of  the  CIPT  phase  space  could  be  accessed  from  the 
far-field.  The  considered  levels  of  nonlinearities  are  available  with 
semiconductor  materials  such  as  GalnAsP  at  frequencies  close  to  the 
band  edge29,  liquid  crystals30  or  multiple  quantum  wells31.  Further¬ 
more,  the  strong  modification  of  the  phase  space  and  the  breaking 
up  of  the  original  invariant  curves  into  islands  is  prerequisite 
behaviour  for  the  onset  of  chaos.  Indeed,  at  very  high  nonlinearities, 
an  ever-increasing  breaking-up  of  the  regular  motion  reaches  a 
point  where  the  phase  space  acquires  strongly  chaotic  characteris¬ 
tics,  as  shown  in  Fig.  3f.  Such  complex  motion  is  directly  related  to 
the  three-body  problem  in  celestial  mechanics  and  underlines  the 
optical-mechanical  analogy  and  its  broad  applicability. 

Methods 

The  invariance  of  Maxwell’s  equations  under  coordinate  transformations  brings 
the  equivalence  between  curved  spacetime  and  local  optical  response  through 
spatially  dependent  permeability  and  permittivity  tensors.  Specifically,  light  motion 
in  the  spacetime  metric  with  the  line  element 

ds2  =  g00dt2  - gu  dx'dx1  (4) 

is  equivalent  to  that  in  an  effective  media  with  local  permeability  and  permittivity 
tensors  given  as 


where  hi  =  Jgii  are  the  Lame  coefficients  of  the  transformation  between  the  curved 
metric  equation  (4)  and  flat  space4,11,  and  the  factor  ^Jg^  is  the  redshift/blueshift 
correction  due  to  time  dilation.  To  eliminate  the  need  for  such  complex  behaviour, 
we  concentrate  on  a  class  of  centrally  symmetric  metrics  that  can  be  written  under 
coordinate  transformation  in  an  isotropic  form  ds2  =g00(r,t )  df2  —g(r,  t) dx2. 
Thus,  the  effective  refraction  index  n  is  also  isotropic  and  can  be  achieved  with 
dielectric  materials  (see  Supplementary  Informations  for  details). 

To  schematically  represent  the  PBH  metrics  (see  Fig.  la),  we  consider  only 
lateral  space  stretching  or  compression  by  means  of  a  normalized  Euclidean 
coordinate  w  =  w(r)/r  as  per  ds2  =  d t2  —  grr(r)  dr2  =  d t2  —  dw2,  where  grr  =  n2. 

In  Fig.  lb,c,  we  also  show  photon  motion  corresponding  to  the  first  integral 
of  the  system  (equation  (1)),  which  is  a  family  of  exponential  spiral  trajectories 
r  =  a  x  ( 1  —  e(<Po-<p)a/byi  ?  where  <p0  is  the  angle  of  approach  at  infinity. 

The  DOM  effective  indices  and  volume  fractions  of  the  inclusions  (Figs  2  and  3) 
are  estimated  from  the  effective-medium  approximation32  and  Maxwell-Garnett 
theory33  for  the  dielectric-dialectic  and  metal-dielectric  composites,  respectively. 
Commercial  software  FEMLAB  has  been  used  for  the  electromagnetic 
finite-difference  frequency  domain  (FDFD)  calculations  throughout  the  paper. 
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